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Executive  Summary 

The  Army  Research  Laboratory- Weapons  and  Materials  Research  Directorate 
(ARL-WMRD)  established  a  Cooperative  Agreement  with  the  University  of  Nebraska- 
Lincoln  (UNL)  to  create  a  Center  for  Research  on  Trauma  Mechanics  and  Blast 
Mitigation.  Research  conducted  under  this  Agreement  was  a  collaborative  effort,  with  the 
ARL-WMRD  and  UNL  sharing  both  leadership  and  research  responsibilities.  Research 
within  the  Center  was  focused  on  development  of  new  materials  and  technologies 
relevant  to  trauma  mechanics,  blast  mitigation  and  weapons  detection.  Specific  studies 
within  the  Center  were  collaborative  in  nature,  with  scientific  personnel  from  ARL 
interacting  closely  with  faculty  at  UNL. 

This  report  summarizes  the  final  results  achieved  for  specific  projects  performed 
under  the  Agreement.  These  projects  had  the  common  objective  of  providing  new 
materials  and  technologies  in  the  areas  of  trauma  mechanics.  The  overall  goal  of  these 
projects  was  to  understand  how  blast  waves  affect  the  head  and  brain  of  a  fully  equipped 
soldier  in  the  field  by  studying  deformation  and  damage.  This  research  also  pursed  to 
couple  macroscopic  effects  with  changes  at  the  microscopic  level  in  cells,  tissues,  and 
organs. 

To  achieve  results  UNL  conducted  a  synergistic  combination  of  experiments, 
modeling,  and  simulation  which  are  unique  approaches.  The  focus  of  the  project  was  on 
wave  propagation  effects  on  the  skull  and  brain  especially  under  mild  TBI  pressure 
loading  conditions.  The  specific  objectives  included: 

1.  Develop  a  state-of-the-art  laboratory  facility  for  generating  well  controlled  air 
shockwaves  that  can  mimic  blast  waves  induced  by  free-field  explosions  and  for 
conducting  well  defined  and  accurately  measured  blast  wave  experiments. 

2.  Study  the  effect  of  blast  waves  or  pressure  pulses  generated  on  a  human  head  with 
and  without  a  protective  helmet  (experiments  and  simulation). 
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3.  Develop  a  multiscale  constitutive  model  of  skull  bone  based  on  experiments  and 
modeling.  Simulate  the  effect  of  pressure  and  impact  loading  on  deformation  and 
intercranial  pressure. 

4.  Investigate  the  functional  and  mechanical  responses  of  neuronal  cells  to  impulsive 
pressure  and  strain. 

This  research  study  was  a  combination  of  experiments  and  modeling  focused  on 
the  biomechanical  response  of  soft  tissue.  The  experiments  were  designed  to  provide 
information  for  the  modeling  efforts  at  a  variety  of  length  scales  and  loading 
configurations  such  that  robust,  comprehensive  models  can  be  developed. 

The  results  of  these  experiments  were  integrated  with  a  multiscale  numerical 
model  developed  at  UNL  and  Sandia  Laboratories.  This  model  simulates  blast  waves  at 
various  distances  and  computes  the  pressure  pulses  experienced  by  the  skull,  brain,  and 
body  extremities. 
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Projects  Conducted 

This  research  program  was  comprised  of  7  collaborative  research  projects.  In  this 
section,  each  of  the  projects  is  addressed  individually  with  an  overall  description  of  the 
project:  objective,  approach,  significance,  and  major  accomplishments. 

A.  Laboratory  Facility  for  Blast  Wave  Simulation  and  Experiments  (Dr.  R.  Feng) 

I.  Research  Objectives 

The  goal  of  this  project  is  to  develop  a  blast  testing  facility  that  can 
reliably  generate  air  shockwaves  with  wave  profiles  mimicking  those  of 
explosion-induced  blast  waves  in  the  mid-  to  far-field  range,  that  ensures  the  test 
sample  (head  form  or  animal  model)  to  be  subjected  to  well  controlled,  accurately 
monitored,  one-dimensional  (ID)  blast  loading,  and  that  produces  time-resolved 
measurements  of  the  test  sample  motion,  deformation  and  pressure  histories  under 
the  blast  loading. 

II.  Approach 

An  explosion-induced  blast  wave  in  open  field  has  a  supersonic  shock 
front  followed  immediately  by  blast  wind  with  time-dependent  pressure  profile 
decaying  from  peak  overpressure  to  peak  underpressure  (vacuum)  and  then 
recovering  back  to  ambient.  Such  a  wave  profile  is  commonly  described  as  the 
Friedlander  wave  profile.  Shock  tubes  have  been  used  since  early  50s  to  generate 
ID  air  shockwaves.  In  this  project,  the  shock  tube  technique  was  adopted  and 
specialized  to  generate  Friedlander  type  of  air  shockwaves.  We  used  pressurized 
nitrogen  or  helium  as  the  driver  gas  and  diaphragms  of  various  designs  to  separate 
the  driver  section  and  the  driven  section  containing  ambient  air.  The  burst 
pressure  of  each  type  of  diaphragm  was  predetermined  experimentally.  The 
sudden  burst  of  a  diaphragm  at  its  designed  burst  pressure  allowed  a  rapid 
expansion  of  the  driver  gas  into  the  driven  section  thus  initiating  a  plane  air 
shockwave  propagating  along  the  shock  tube.  Unique  designs  of  square  tube  and 
adjustable  end  reflector  were  developed  and  implemented  to  enable  the  use  of 
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high-speed  digital  cameras  for  measurement  and  imaging  and  to  ensure  the 
desired  (Friedlander)  wave  loading.  For  diagnostics  and  control,  the  pressure 
profiles  of  the  shockwaves  at  various  locations  along  the  shock  tube  were 
measured  in  time-resolved  manner  with  dynamic  pressure  sensors. 

III.  Significance 

The  blast  testing  facility  is  a  critical  need  for  blast-induced  traumatic  brain 
injury  (bTBI)  research  in  that  it  enables  experimental  studies  of  the  physics  and 
mechanics  underlying  the  transformation  of  blast  wave  incidence  to  brain  strain 
and  pressure,  measurements  of  the  dynamic  response  of  skull  and  brain  tissues 
and  brain  cells  under  blast  loading  condition,  animal  model  experiments  to 
determine  the  blast  intensity  thresholds  and  pathological  mechanisms  for  bTBI  in 
animal  samples,  and  experimental  evaluations  of  the  performances  of  head 
protection  systems  and  designs  in  the  laboratory  environment. 

IV.  Accomplishments 

In  this  period,  we  designed  and  constructed  two  shock  tubes:  a  4” 
cylindrical  uniform  shock  tube  with  variable  driver  and  driven  section  lengths  for 
pulse  shaping  [Fig.  1(a)]  and  a  9”  square  shock  tube  with  optical  windows  for 
digital  imaging  [Fig.  1(b)].  A  28-inch  square  shock  tube  with  variable  driver  and 
driven  section  lengths  and  test  section  position  for  wave  profile  control  was  also 
designed.  A  data  acquisition  system  for  up  to  24  channels  of  synchronized  time- 
resolved  measurements  at  2.5  MHz  frequency  and  the  computer  program  for  its 
remote  control  were  developed.  It  allows  pressure  profile  measurements  at 
multiple  locations  to  accurately  track  the  evolution  of  blast  wave  profile.  A  gas 
handling  system  for  driver  pressurization  and  breech  firing  and  the  network  and 
computer  program  for  its  remote  monitoring  and  control  was  developed. 
Capability  for  generating  well  controlled  highly  repeatable,  shockwaves  with  the 
Friedlander  wave  profile  was  established  for  both  the  4”  and  9”  shock  tubes. 


5 


Fig.  1  :  Blast  wave  simulation  and  testing  facility,  (a)  4”  cylindrical  shock  tube;  (a)  9” 
square  shock  tube;  (c)  9”  tube  transition  section;  (d)  Pressure  sensor  array;  (e)  9”  tube 
adjustable  end  reflector;  and  (f)  ARAMIS  three-dimensional  digital  imaging  system. 

It  was  found  from  the  results  of  the  Experiments  with  the  4”  cylindrical 
uniform  shock  tube  that  the  independently  measured  blast  wave  front  speed  and 
peak  pressure  correlate  very  well  with  ID  shock  wave  theory  for  ideal  gas  giving 
proof  that  the  simulated  blast  wave  is  a  ID  plane  wave.  The  device  can  be  used 
for  sensor  calibration. 

The  unique  9”  square  shock  tube  has  a  conical  transition  section 
expanding  the  cross  section  area  from  that  of  a  4”  cylindrical  breech  to  that  of  a 
9”  square  tube  [Fig.  1(c)].  Although  the  air  flow  in  and  near  the  transition  section 
is  complex  and  multi-dimensional,  it  was  validated  experimentally  that  the 
resulting  shockwave  in  the  test  section  (8  ft  away  from  the  transition  section)  is 
essentially  a  ID  blast  wave  in  the  middle.  By  incorporating  an  adjustable  end 
reflector  at  the  shock  tube  exit,  we  found  from  experimental  measurements  that 
the  amplitude  and  duration  of  the  underpressure  part  of  simulated  blast  wave  can 
be  well  controlled  with  a  right  combination  of  the  breech  length  and  the  end 
reflector  position.  A  comparison  of  the  pressure  wave  profiles  with  and  without 
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Fig.  2:  Comparison  of  pressure  wave  profiles  with  and  without  adjustable  end  reflector. 

the  adjustable  end  reflector  is  shown  in  Fig.  2.  We  also  successfully  incorporated 
a  state-of-the-art  ARAMIS  high-speed  three-dimensional  (3D)  digital  imaging 
correlation  system  into  the  blast  testing  with  the  9-inch  shock  tube  [Fig.  1(f)].  The 
synchronized  digital  imaging  and  3D  imaging  correlation  technique  can  provide 
visualization  of  the  test  sample  motion  under  the  blast  loading  and  more 
importantly  time-resolved  measurements  of  the  sample  displacement  and  surface 
strain  evolutions  and  therefore  add  a  unique  capability  to  our  blast  simulation  and 
test  facility. 

With  the  9”  shock  tube,  we  carried  out  a  series  of  blast  tests  on  tubular 
cylindrical  polycarbonate  samples  of  various  sizes  and  with  various  fills  (water, 
mineral  oil  and  silicone  gel).  In  these  experiments,  we  used  a  high-precision 
piezoelectric  sensor  array  to  monitor  the  blast  wave  profile  before  and  after  the 
interaction  with  the  test  sample,  the  ARAMIS  3D  digital  imaging  correlation 
system  to  measure  the  sample  displacement  and  surface  strain  histories  under  the 
blast  loading,  and  FISO  fiber  optics  sensor  to  measure  the  resulting  pressure 
history  in  the  material  filled  the  sample.  Presented  in  Fig.  3  are  the  experimental 
set-up  and  a  set  of  ARAMIS  measurements. 
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Fig.  3:  Blast  wave  experiment  on  fluid-filled  polycarbonate  tube,  (a)  Experimental  set-up; 
(b)  displacement  history  at  point  0  (near  of  tube  side  surface  center);  (c)  Vertical  strain 
history  at  point  0;  (d)  Snap  shot  of  vertical  strain  field  at  the  time  marked  with  red  line  in 
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B.  Realistic  Explosive  Dummy  (RED)  Head  Development  (Dr.  Nelson) 

I.  Research  Objectives 

The  main  objective  of  this  project  is  to  create  a  surrogate  headform  whose 
mechanical  behavior  is  similar  to  that  of  the  human  head.  Through  shock-tube 
testing  of  the  headform  and  cross-validation  using  computational  simulations,  we 
can  learn  more  about  the  possible  mechanisms  of  blast-induced  traumatic  brain 
injury  and  the  high-speed  interactions  that  occur  inside  the  head  during  such  an 
event. 

II.  Approach 

Most  existing  headforms  are  designed  to  replicate  the  amount  and 
distribution  of  mass  in  the  head  but  not  to  replicate  the  actual  mechanical 
behavior  of  tissues  in  the  head. 

More  advanced  headforms  may  give  special  attention  to  specific 
structures,  such  as  the  eyes,  and  include  some  force  or  pressure  sensors. 

However,  no  headform  had  previously  included  any  truly  brain-like  structure  or 
means  of  measuring  high-speed  pressure  changes  inside  the  skull  cavity.  Through 
study  of  the  mechanical  properties  of  the  various  tissue  components  in  the  human 
head,  we  were  able  to  create  a  headform  with  mechanical  behavior  duplicating  to 
a  significant  degree  that  of  the  human  head,  including  the  brain.  We  also 
implemented  a  high-speed  data  capture  system  with  various  types  of  sensors  on/in 
the  headform  in  order  to  record  the  phenomena  occurring  during  shock  or  blast 
events. 

III.  Significance 

The  main  significant  contribution  of  this  work  is  the  ability  to  study  the 
high-speed  phenomena  in  the  human  head  which  may  contribute  to  the 
development  of  traumatic  brain  injury  in  individuals  exposed  to  blasts.  This  has 
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direct  relevance  to  warfighter  protection,  care,  diagnosis,  and  so  forth.  It  also  has 
similar  civilian  applications  for  industrial  accidents. 

By  gaining  a  more  in-depth  knowledge  about  these  phenomena,  one 
important  potential  direction  is  development  of  improved  protective  equipment. 

IV.  Accomplishments 

Two  iterations  of  a  biofidelic  headform  (RED  Head)  have  been  developed 
with  high-speed  instrumentation  to  measure  pressure,  strain,  etc.  Shock-tube  tests 
have  been  performed  and  data  collected  and  analyzed.  Ongoing  work  includes 
further  characterization  of  the  effects  of  material  and  geometry  on  shock  wave 
propagation  through  these  types  of  shock  targets,  including  additional 
experimental  testing. 
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C.  TBI-Materials  Modeling  and  Characterization  and  3-D  Optical  Measurements 
(Dr.  Negahban) 


I.  Research  Objectives 

1)  Modeling  and  Characterization 

The  materials  modeling  and  characterization  group  has  been  working  on 
characterizing  the  skin,  skull  and  brain,  and  the  associated  natural  and  man-made 
materials  used  in  the  design  and  construction  of  the  RED-Head,  in  protective 
components,  and  the  blast  facility.  The  experimental  setup  and  numerical  analysis 
within  the  TBI  blast  community  and  the  UNL  Blast  facility  need  accurate  material 
properties  and  models  for  the  different  components  that  compose  the  human  and 
animal  testing  setups,  the  dummy  RED-Head,  and  protective  armor.  The  group 
works  to  collect  data,  develop  new  experimental  techniques,  and  construct  models 
to  better  support  the  communities  ability  simulate,  design,  and  develop  systems  to 
protect  against  TBI. 

2)  3-D  Optical  Measurements 

The  3-D  optical  measurement  group  works  with  provide  continuous 
surface  strain  measurements  capabilities  and  develop  this  technique  for  dynamic 
testing  system.  In  doing  so  the  group  supports  the  RED-Head  group,  the  dynamic 
testing  group. 

II.  Approach 

Since  the  modeling  and  characterization  group  worked  to  support  many 
different  constituents,  the  approach  taken  included: 

(a)  In  support  of  the  RED-Head  development  and  materials 
characterization  we  conducted  dynamic  mechanical  testing  (DMA)  on  materials 
being  considered  and  used  in  the  RED-Head.  This  included  DAM  results  in 
compression,  double  shear,  and/or  cantilevered  beam  modes  for  the  polymer  skin 
and  skull  materials  provided  by  the  manufacturer  of  the  RED-Head  and  for 
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different  silicone  gels  and  other  materials  considered  for  the  brain.  Similar  tests 
were  conducted  on  pig  brain,  liver,  and  cow  hip  bone  to  provide  a  comparative 
measure  between  the  artificial  materials  used  in  the  RED-Head  and  biological 
materials  (Fig.  4);  (b)  In  an  effort  to  consider  how  to  improve  the  skull  model  to 
better  mimic  the  actual  skull,  we  did  an  experimental  study  of  the  effects  of 
different  microstructures  on  the  average  ultrasonic  signal  transmitted  through 
polymer  samples  with  different  patterns  of  voids  constructed  using  rapid 
prototyping;  (c)  In  an  effort  to  better  characterize  the  bone  structure  in  the  skull, 
we  conducted  a  study  of  the  anisotropic  properties  of  bovine  hip  bone  taken  along 
different  directions  using  both  DMA  and  ultrasonic  waves;  (d)  In  an  effort  to 
overcome  the  strong  dissipation  that  controls  the  possibility  of  measurements  in 
soft  tissue,  and  in  particular  the  brain,  we  have  been  developing  a  method  to 
characterize  the  wave  speed  and  attenuation  under  load  using  reflected  signal 
changes  at  different  frequencies;  (e)  To  improve  models  currently  used  in  the 
simulation  of  the  skin,  skull  and  brain,  we  constructed  viscoelastic  models  to 
characterize  both  the  shear  and  dilatational  response  of  these  components  under 
infinitesimal  deformations  using  Prony  series  models  of  the  relaxation;  (f)  To 
better  capture  the  dynamic  response  of  the  brain  we  developed  a  method  to  fit  a 
continuous  relaxation  function  to  shear  data,  allowing  a  much  better  fit  of  both  the 
storage  and  loss  moduli;  (g)  To  capture  the  nonlinear  response  of  the  skin,  skull 
and  brain  under  the  more  realistic  nonlinear  setting,  we  have  been  developing  a 
theory  for  wave  propagation  in  nonlinear  viscoelastic  materials. 

Since  the  3-D  optical  measurement  group  work  in  both  development  of 
technique  and  support  of  current  activities  by  other  groups  working  on  the  TBI 
project,  the  work  took  different  approaches  to  the  different  responsibilities  that 
included  working  with  the  RED-Head  group,  the  dynamic  testing  group,  and  the 
blast  facility  to  conduct  low  and  high  speed  optical  measurements  in  support  of 
the  following  activities:  (a)  Support  of  the  RED-Head  in  conducting  dynamic 
testing  on  different  gels  for  use  as  the  brain  simulation  material;  (b)  Dynamic 
characterization  of  the  blast  facility,  including  the  characterization  of  the  effect  of 
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windows,  characterizing  the  shock  induced  change  in  the  index  of  diffraction  of 
the  gas,  and  general  high  speed  measurements  during  development  of  the  facility; 
(c)  Dynamic  characterization  of  sample  surface  deformation  inside  and  outside 
the  blast  tube  for  both  the  RED-Head  and  simple  sample  geometries  used  to 
calibrate  the  gun,  the  measurement  systems  and  the  computer  simulations;  (d) 
High  speed  characterization  on  the  Kolsky  bar  and  on  other  setups  (Figs.  5  and  6). 


Selection  of  materials  for  the  RED  head 
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Fig.4:  Comparing  model  materials  and  true  response. 


III.  Significance 

The  constitutive  characterization  and  modeling  group  both  established 
support  the  design,  development,  and  analysis  TBI  facility  and  worked  on  the 
development  of  new  characterization  methods  for  soft  tissues  and  large 
deformations. 


The  3-D  optical  measurement  group  developed  optical  characterization 
competencies  for  quasi-static  to  dynamic  surface  strain  measurements  using  an 
ARAMIS  low  and  high  speed  system,  and  supported  the  design  and  the 
development  of  the  blast  facility. 
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Fig.  5:  3-D  Optical  measurement  for  characterization  of  blast  wave  distortion. 


14 


Fig.  6:  3-D  Optical  measurement  for  characterization  of  blast  interaction  with  a  cylinder. 

IV.  Accomplishments 

(1)  The  group  established  mechanical  characterization  between  human  and 
animal  head  components  including  skin,  bone,  and  brain,  and  their 
counterparts  used  in  the  RED-Head. 

(2)  The  group  established  improved  modeling  for  the  characterization  of  skin, 
skull,  and  brain  viscoelastic  properties. 

(3)  The  group  improved  the  modeling  of  the  linear  viscoelastic  modeling  of 
head  components  by  developing  a  continuous  response  spectrum  (Fig.  7). 

(4)  The  group  worked  on  characterization  techniques  for  obtaining  elastic  and 
viscoelastic  properties  of  highly  attenuating  materials  such  as  biological 
tissue  using  wave  reflection  methods  and  conducted  preliminary 
evaluations. 
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(5)  The  group  worked  on  perturbation  in  nonlinear  viscoelastic  loaded 
materials  to  establish  a  better  characterization/modeling  method  for 
biological  materials  undergoing  large  deformations. 


Modeling  of  Soft  Materials 


Improvement  of  current  models  by  using  extended  Prony  series  and  continuous 
relaxation  models: 


o®  =  n  crs  (t)  ds 

Continuous  element 
model  provides 
better  control  over 
data  fitting  without 
difficulty  in  fitting  both 
relaxation  and 
storage  modulus 


Circular  Frequency  (rad/s) 


Data  from  Nicolle  et.  all  (2005) 


Lincoln 


Fig.  7:  Improving  modeling  by  using  continuous  spectra. 
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D.  Influence  of  Skull  Microstructure  on  the  Profile  of  Pressure  Waves  (Dr.  Turner) 


I.  Research  Objectives 

The  objective  of  this  project  is  to  quantify  the  influence  of  skull 
microstructure  on  the  profile  of  pressure  waves  that  pass  through  the  skull. 

II.  Significance 

The  flat  bones  of  the  skull  are  comprised  of  the  outer  and  inner  tables 
(compact  bone)  and  the  diploe  (porous  bone)  on  the  interior.  This  variation  in 
microstructure  within  the  skull  causes  a  disruption  in  the  coherent  nature  of  waves 
that  are  incident  on  the  skull.  Blasts  that  result  from  IEDs  and  other  small 
ordnances  are  often  characterized  by  a  fast  rise  time  such  that  they  are  fairly  broad 
band  (up  to  ~1  MHz  in  frequency  content). 

The  scattering  behavior  that  will  result  from  the  skull  microstructure  will 
reduce  the  pressure  profile  (i.e.,  strain)  and  alter  the  frequency  content  (i.e.,  strain 
rate)  of  the  incident  wave.  Currently,  very  little  is  known  about  the  cells  of  the 
brain  with  respect  to  their  response  to  excitations  that  include  a  wide  frequency 
range  since  strain  rate  effects  on  single  cells  or  groups  of  cells  remains  largely 
unexplored.  Thus,  a  clear  understanding  of  the  influence  of  skull  microstructure 
on  wave  transmission  is  needed.  As  part  of  the  proposed  research,  a  variety  of 
techniques  will  be  used  to  quantify  this  behavior  and  connect  the  skull 
architecture  to  the  wave  behavior.  Human  skull  samples  will  be  used  as  the  basis 
for  all  measurements  so  that  the  results  will  be  directly  applicable  to  the 
warfighter. 

III.  Accomplishments 

Previous  research  on  this  aspect  of  the  TBI  research  was  started  in  the  last 
two  years  at  UNL  in  the  Pi’s  group.  A  human  skull  sample  was  obtained  from 
UNMC  in  order  to  begin  exploration  of  this  research  area  that  has  been  largely 


17 


neglected.  The  sample  was  obtained  from  a  healthy  78-year-old  male  patient  who 
required  brain  surgery  (Fig.  8). 

The  skull  piece  could  no  longer  be  replaced  after  the  surgery  and  was 
provided  for  this  research.  Small  samples  were  cut  from  the  larger  sample  so  that 
microstructure  and  property  information  could  be  obtained.  MicroCT  and 
nanoindentation  measurements  were  made  on  the  four  samples  to  quantify  these 
aspects  of  the  skull  response.  Porosity  ranged  from  17-30  %  and  the  storage 
modulus  was  in  the  range  of  20  GPa.  The  optical  images  were  used  to  construct 
2D  finite  element  (FE)  models  so  that  an  initial  understanding  of  the  influence  of 
microstructure  could  be  obtained. 

The  FE  models  were  developed  using  ABAQUS  software  and  are 
highlighted  in  Fig.  2.  The  model  consisted  of  a  region  that  represented  the  skull 
(either  as  solid  or  porous)  and  a  region  that  represented  the  brain  (modeled  as  a 
fluid).  Different  pressure  profiles,  exciting  the  entire  top  surface  of  the  model, 
were  explored  as  incident  wave  inputs  including  a  three-cycle  Gaussian  pulse  and 
a  profile  representative  of  a  blast  pulse. 
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Fig.  8:  Example  skull  sample  obtained  from  UNMC  (Dr.  P.  Lennarson).  Four  samples 
were  cut  from  the  larger  sample.  The  smaller  samples  were  used  to  get  microstructure 
information  (using  microCT)  and  viscoelastic  property  information  (using 
nanoindentation) . 
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Top  surface  excited  with  Gaussian  pulse 

Skull  (solid) 


Brain  (fluid) 

Average  pressure  profile  alongline  within 
brain  shows  multiple  reflections  from  free 
surfaces  when  skull  is  homogeneous. 


When  the  skull  is  heterogeneous  the 
pressure  profile  in  the  brain  is  greatly 
reduced  and  clearly  shows  scatteringdue 


to  porosity. 


Line  average 


Fig.  9:  Example  finite  element  results  that  highlight  the  potential  impact  skull 
microstructure  can  have  on  the  transmission  of  pressure  through  the  skull. 
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E.  Fluid-Structure  Interaction  (FSI)  (Dr.  Bobaru) 


I.  Research  Objective 

Many  of  Traumatic  Brain  Injury  (TBI)  phenomena  involve  complex  Fluid- 
Structure  Interaction  (FSI)  induced  by  blast  waves  propagation.  The  objective  of 
this  part  of  the  project  was  to  find  and  develop  numerical  verification  methods  for 
the  simulation  of  the  FSI  problem  to  better  understand  the  mechanical  behavior  of 
a  three-dimensional  test  sample  in  a  shock  tube  experimental  setup.  The  goal  was 
to  model,  at  the  same  scale  and  geometry  as  the  shock  tube  developed  and  used  in 
experiments  performed  by  another  group  (Dr.  Feng)  on  the  same  project.  The 
computational  model  we  provided  helped  the  experimental  group  to  design  more 
reliable  experiments  and  improve  the  actual  shock  tube.  For  example,  better 
location  of  the  sensors  on  the  sample  tested  was  found  using  our  computational 
models. 

II.  Approach 

The  main  target  model  was  the  UNL  Shock  Tube  facility  (Dr.  Feng 
group).  We  modeled  the  propagation  of  blast  waves  using  the  same  geometry  and 
scale  as  the  actual  shock  tube  and  using  the  same  blast  wave  generation 
mechanism  that  induced  the  blast  wave  in  experiments.  We  first  verified  the  main 
characteristics  of  blast  wave  propagations,  such  as  the  shape  of  blast  wave  profile, 
peak  pressure  magnitude,  etc.  We  employed  the  finite  volume  method  built  in 
ABAQUS.  We  constructed  the  model  for  the  blast  wave  generation  mechanism 
and  implemented  this  into  ABAQUS.  Next,  we  modeled  the  coupled  FSI  problem 
by  coupling  the  immersed  finite  element  method  (for  the  solid  materials)  with  the 
finite  volume  method  (for  the  fluid  part).  Many  challenging  issues  around  the 
intermediate  region  shared  by  both  solid  and  fluid  elements  arose:  inconsistent 
boundary,  contact  condition,  etc.  We  found  solutions  to  these  issues.  We 
investigated  the  correlations  of  the  fluid  pressure  magnitude  and  profile  ahead  and 
around  the  solid  obstacle  and  compared  with  experiments. 
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III.  Significance 

Understanding  the  fluid-structure  interaction  in  blast  conditions  is 
paramount  to  progress  in  the  TBI  problems.  The  design  of  experiments  in  the 
shock  tube  used  the  results  we  produced  for  their  calibration  and  testing. 

IV.  Accomplishments 

In  this  small-scale  portion  of  the  project  we  have  developed  the 
generalized  Fluid-Structure  Interaction  (FSI)  computational  model  in  ABAQUS 
and  performed  computations  for  the  experiments  conducted  by  Dr.  Feng’s  group 
at  UNL’s  shock  tube  facility. 
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F.  Effects  of  Impulsive  Pressurization  on  Brain  Cell  Functions  (Dr.  Lim,  Dr.  Feng) 


I.  Research  Objectives 

A  blast-induced  traumatic  brain  injury  (bTBI)  caused  by  exposure  to  blast 
of  various  explosive  devices,  including  improvised  explosive  device  (IED),  may 
result  in  post-concussion  symptoms  that  are  easily  overlooked  immediately  after 
the  event  but  pose  long  term  health  problems.  A  better  understanding  of  the 
cellular  level  mechanisms,  through  which  an  IED  blast  induces  the  primary  bTBI 
and  its  secondary  progression,  is  critically  important  for  developing  effective 
protection  gears,  early  diagnostics  techniques/procedures,  and  proper  treatments 
against  bTBI. 

II.  Approach 

To  investigate  the  role  of  dynamic  short-duration  overpressure  (one  of  the 
speculated  factors)  in  the  brain  cell  response  to  bTBI  relevant  blast  loading,  we 
have  developed  novel  impulsive  brain  cell  pressurization  equipment  based  on 
Kolsky  bar  set-up.  This  apparatus  allows  examination  of  brain  cell  functions  after 
exposure  to  in  vitro  conditions  mimicking  blast-induced  pressure  loading. 

III.  Significance 

Currently,  bTBI  is  measured  using  a  combination  of  magnetic  resonance 
imaging  (MRI)  and  the  assessment  of  the  duration  of  loss  of  consciousness  or 
amnesia.  However,  no  technology  currently  exists  for  reliably  replicating  the  blast 
loading  of  the  brain  cells  in  vitro.  While  the  mechanism  of  how  blast  waves  affect 
the  brain  and  brain  cells  is  not  fully  revealed,  collective  data  suggest  that  bTBI  is 
a  short-duration  overpressure-induced,  stress  wave-dominant  phenomenon. 
Therefore,  the  development  of  technology  that  enables  the  evaluation  of 
impulsive  overpressure  effects  on  brain  cells,  in  both  short  and  long-term 
perspectives,  under  bTBI  relevant  blast  loading  conditions  will  provide  a  huge 
opportunity  to  reveal  the  nature  of  bTBI. 
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IV.  Accomplishments 

A  technique  to  generate  fast-rising,  long-duration,  single  compressive 
stress  pulses  has  been  developed  by  incorporating  a  new  axial  loading  system 
(which  consists  of  a  two-way  hydraulic  piston  and  an  inertia  stopper)  in  a  Kolsky 
torsion  bar  device  available  in  Feng’s  lab  (Laboratory  for  Dynamic  Materials 
Characterization).  A  cell  pressurization  chamber  that  can  convert  axial  stress 
pulse  imposed  on  its  exterior  to  well  defined  pressure  pulse  in  the  culture  medium 
fluid  contained  in  its  interior  has  been  developed.  The  fixation  of  cell  cultured 
glass  slide  within  the  chamber  has  been  accomplished.  The  top  view  of  Fig.  10 
shows  the  Kolsky  bar  cell  pressurization  set-up  and  schematic  of  force  transfer 
during  the  pressurization.  Using  this  apparatus,  single  compressive  stress  pulses 
with  a  rise  time  of  approximately  0.05  ms,  a  pulse  duration  of  approximately  0.7 
ms,  and  stress  levels  at  0.5-2  MPa  have  been  accomplished  as  a  transient 
overpressure  stimuli  to  the  cultured  neuronal  cells.  The  bottom  view  of  Fig.  10 
shows  an  example  that  0.5,  1,  and  2  MPa  level  impulsive  pressurizations. 

Human  SH-SY5Y  human  neuroblastoma  cell  (ATCC,  CRL-2266)  culture 
and  pressurization  protocols  have  been  established.  Before  pressurization,  SH- 
SY5Y  cells  were  treated  with  10  pM  retinoic  acid  (RA)  for  7  days  for 
neurite/axonal  outgrowth.  Cells  treated  with  RA  the  same,  placed  inside  the 
pressurization  vessel  the  same,  but  not  pressurized  served  as  an  unpressurized 
chamber  control.  Effects  of  bTBI  relevant  impulsive  pressurization  on  the 
cytoskeletal  structure  of  neuronal  cells  were  examined  (Fig.  11).  Two  of  the 
important  cytoskeletal  structures,  actin  (red)  and  tubulin  (green),  were  clearly 
observed  in  the  unpressurized  control  cells  following  the  neurite/axonal  structure. 
Neuronal  cells  exposed  to  impulsive  pressurization  at  2  MPa  magnitudes  and  0.7 
ms  duration  showed  distinct  breakdown  of  the  neurite/axonal  structure,  as  seen  by 
shorter  or  no  actin/microtubule  staining.  Microtubules  serve  not  only  as 
cytoskeletal  components  but  also  play  a  role  as  a  marker  of  neurogenesis. 
Microtubles  are  heavily  involved  in  the  expression  of  microtubule-associated 
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protein-2  (MAP2),  which  is  a  distinct  marker  of  neuronal  differentiation.  Our  data 
thus  suggest  that  bTBI-relevant  pressurization,  that  induces  changes  in  apparent 
neurite  morphology  and  cytoskeletal  components,  may  further  affect  neuronal  cell 
function  and  activity. 


Fig.  10:  Schematic  of  impulsive  cell  pressurization  device  based  on  Kolsky  bar. 
The  bottom  diagram  shows  schematic  of  force  transfer  during  the  pressurization. 

We  could  control  the  impulsive  pressurization  levels  at  0.5-2  MPa 
magnitude.  The  pressures  recorded  with  a  Fiso  pressure  gauge  had  the  same 
square  pulse  shape  as  the  transmitted  bar  measurements. 
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Fig.  11:  Human  neuronal  cells  pressurized  at  2  MPa  and  0.7  ms  duration  display  distinct 
neurite/axonal  breakdown  relative  to  chamber  control  cells,  e.g.,  shortened  or  removed 
neurite/axons  in  the  immunofluorescence  of  actin  (red)  and  a  tubulin  (green). 


26 


G.  Computational  Modeling  of  Human  Head  and  Testing  of  Neurons  (Dr.  Chandra) 


I.  Research  Objectives 

1)  Modeling,  simulation  of  the  head  (brain)  of  war  fighters  (humans)  with 
high  biofidelic  anatomical  details. 

2)  Examine  the  role  of  protective  helmets  in  mitigating  the  effect  of  shock- 
blast  induced  by  EEDs. 

3)  Mechanisms  of  deformation  and  damage  of  in-vitro  neurons  and  high  rates 
of  loading. 

II.  Approach 

Modeling  and  simulations:  Integrate  a  series  of  medical  MRI  and  CT 
images  and  construct  a  3D  brain-head  model.  The  model  is  anatomically  accurate; 
and  can  be  discretized  into  finite  elements  for  analysis  and  validation  with 
experimental  results. 

Role  of  protective  helmets:  Digitize  accurately  helmets  of  various 
generations  and  origin  and  import  into  numerical  models.  Include  the  effects  of 
harness  and  various  pads  with  high  geometric  accuracy.  The  effect  of  helmets  can 
be  analyzed  by  understanding  the  mechanics  of  shock-blast  interactions  with  the 
head. 

In-vitro  neuron  testing:  Culture  human  neurons  in  the  laboratory  and  test 
them  under  different  strain-rate  and  establish  strain-response  curves.  The  response 
may  be  biochemical  alterations  of  proteins  or  necrosis. 

III.  Significance 

Blast  TBI  is  one  of  the  critical  unresolved  issues  that  have  a  political, 
economic  and  personnel  consequences  to  US  army.  Understanding  if  and  why 
blasts  cause  mild  TBI  is  very  timely.  Once  mechanisms  are  identified  then 
mitigative  strategies  can  be  formulated. 
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IV.  Accomplishments 


Significant  progress  has  been  made  in  the  modeling  aspects.  We  now  can 
build  accurate  geometric  and  biofidelic  models  and  simulate  shock-blast 
interactions. 


Technology  Transfer 

I.  The  ‘RED  Head’  or  ‘Realistic  Explosive  Dummy  Head’  development  project 
focuses  on  head  injuries  caused  by  shock  waves  generated  by  explosions.  To 
develop  such  a  head  model,  materials  analysis  was  performed  to  find  synthetic 
materials  to  simulate  human  head  components  such  as  skin,  skull  and  brain.  The 
geometries  of  these  components  were  carefully  controlled  to  simulate  it  as  a 
human  head  surrogate.  The  head  model  was  instrumented  with  various  sensors  to 
understand  the  attenuation  of  the  pressure  waves  through  the  brain  and  the  flexure 
of  the  skull.  The  prototypes  developed  as  part  of  this  project  are  available  for  loan 
to  and  use  by  DOD  or  other  personnel  if  desired.  The  details  of  their  design  can 
also  be  shared  and  have  been  presented  to  some  extent  in  publications,  and 
presentations  which  have  been  shared  in  the  2009  summer  progress  meeting  with 
DOD  personnel. 

II.  The  numerical  models  developed  in  ABAQUS  can  be  used  by  researchers  at  DoD 
to  better  understand  the  complex  behavior  of  structures  and  systems  under  blast 
conditions. 

III.  Throughout  the  funding  period,  constant  interactions  with  US  army  (ARO,  PEO, 
ARL,  US  Army  Medical  command)  occurred.  The  PI  Dr.  Namas  Chandra  is  a  part 
of  an  expert  panel  organized  by  Michael  Leggiri.  He  also  presented  his  findings  in 
the  army  meeting  at  St.  Petersburg  in  June  2010. 
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Abstract 

Blast  induced  traumatic  brain  injury  (bTBI)  is  signature  injury  in  recent  combat 
scenarios  involving  improvised  explosive  devices  (IEDs).  The  exact  mechanisms  of  bTBI 
are  still  unclear  and  protective  role  of  helmet  and  body  armor  is  often  questioned. 

Because  of  the  complex  and  highly  nonlinear  fluid-structure  interactions  involved  in 
bTBI  and  the  limitation  of  diagnostic  methods  in  comparison,  robust  computational 
modeling  with  experimental  validation  is  necessary  in  order  to  understand  mechanisms  of 
bTBI.  In  an  effort  to  understand  the  mechanics  of  blast-induced  traumatic  brain  injury 
(bTBI),  we  have  built  numerical  head  models  utilizing  the  arbitrary  Euler-Lagrangian 
coupling  method  available  in  the  ABAQUS  finite  element  code.  Our  simulation 
methodology  allows  accurate  concurrent  simulations  of  the  formation  and  propagation  of 
blast-like  air  shockwave,  the  fluid- structure  interactions  between  the  shockwave  and  the 
head,  and  the  stress  wave  propagation  within  the  head.  In  this  method,  an  Eulerian  mesh 
is  used  to  model  the  air  surrounding  the  head  and  a  Lagrangian  mesh  is  used  to  model  the 
head.  The  shock  wave  is  generated  by  releasing  high  pressure,  high  temperature 
compressed  air  into  the  atmospheric  air  at  the  beginning  of  the  analysis.  The  predefined 
pressure  and  temperature  of  compressed  air  is  adjusted  so  as  to  generate  blast 
overpressures  equivalent  to  blast  overpressures  in  nonlethal  blast  injury,  threshold  lung 
injury  and  50%  lethal  lung  injury.  We  have  built  two  types  of  numerical  head  models. 

The  first  head  model  is  surrogate  dummy  head  model  consisting  of  skin,  skull,  brain  and 
the  neck.  The  material  properties  for  various  head  components  such  as  dummy  skull, 
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dummy  skin  and  brain  simulant  are  obtained  from  various  tests  performed  in  house.  This 
dummy  head  model  is  validated  against  shock  tube  experiments  conducted  in  in-house 
shock  testing  facility.  This  model  has  served  as  'benchmark'  for  our  modeling  effort  and 
suggests  reliability  of  our  modeling  effort.  Second  head  model  is  anatomically  detailed 
three  dimensional  head  model  (real  head  model  from  here  on  )  developed  from 
segmentation  of  MRI  images  available  from  Visible  Human  Project  of  the  National 
Library  of  Medicine  .  This  model  consists  of  anatomical  details  of  head  such  as  skull, 
ventricle,  white  matter,  grey  matter,  CSF,  air  sinus,  tentorial  membranes  and  falx.  The 
material  properties  of  various  tissues  of  the  real  head  model  are  either  taken  from 
literature  or  obtained  from  in-house  testing.  This  real  head  model  is  subjected  to  various 
shock  intensities  under  various  head  orientations.  The  head  response  under  various  shock 
intensities  and  different  head  orientations  is  studied  and  compared  on  the  basis  of 
concussive  and  axonal  injury  criterions  available  in  the  literature.  Parametric  studies  over 
range  of  material  properties  are  conducted  and  critical  components  and  areas  of  head 
under  blast  loadings  are  identified. 


Session:  2-9-1  Modeling  of  Brain  Injury 
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EXTENDED  ABSTRACT 


ABSTRACT 

A  numerical  simulation  (finite  element  method  (FEM))  is 
used  to  determine  the  local  mechanical  response  of  a 
neuron  based  on  the  real  neuron  geometry.  The  real  3D 
geometry  of  a  neuron  can  be  constructed  from  the  2D 
confocal  image  stack  which  is  taken  per  micron  along  the 
thickness  direction  from  the  top  surface  of  neuron  to  the 
substrate.  By  matching  the  simulated  deformation  of  the 
neuron  with  the  measured  results  from  the  confocal 
microscopy  images,  the  relative  strength  of  the  neuron 
components  can  be  determined  based  on  the  reverse 
analysis  of  FEM.  Then,  the  neuron  local  mechanical 
response  and  the  relationship  between  the  local  response 
and  the  global  applied  deformation  can  be  obtained, 
which  will  be  the  first  step  of  building  the  mechanical 
injury  model  of  neuron. 

INTRODUCTION 

The  mechanism  of  mild  traumatic  brain  injury  (mTBI)  is 
directly  related  to  the  relationship  between  the  mechanical 
response  of  neurons  and  their  biological/chemical 
functions  since  the  neuron  is  the  main  functional 
component  of  brain.[l]  Therefore,  understanding  the 
mechanical  response  of  neurons  is  an  important  first  step 
to  understand  the  mechanism  of  mTBI.  Typically,  the 
mechanical  response  of  neurons  is  investigated  based  on 
the  deformation  of  in  vitro  model,  in  which  the  substrate 
deformation  is  considered  to  be  the  deformation  of 


neurons.  However,  neurons  have  the  irregular  shape  and 
mainly  include  cell  main  body,  neurite  and  axon  which 
have  highly  different  geometries  and  mechanical 
properties.  Under  the  same  global  deformation  (substrate 
deformation,  e.g.  substrate  stretching),  the  different  parts 
of  a  neuron  might  have  the  highly  different  local 
deformations.  Therefore,  it  is  important  to  understand  the 
local  mechanical  response  of  neurons  in  order  to  get  the 
accurate  mechanical  injury  model  of  the  neuron  which 
describes  the  relationship  between  the  neuron  mechanical 
response  and  its  function  loss. 

COMPUTATIONAL  METHODS 

A  numerical  simulation  (finite  element  method  (FEM)) 
can  be  used  to  determine  the  local  mechanical  response  of 
a  neuron  based  on  the  real  neuron  geometry.  In  FEM 
simulations,  the  real  neuron  geometry  can  be  considered 
in  both  2D  and  3D.  With  the  help  of  the  software 
(Mimics), [2]  the  real  3D  geometry  of  a  neuron  can  be 
constructed  from  the  2D  confocal  image  stack  which  is 
taken  per  micron  along  z-direction  from  the  substrate  to 
the  top  surface  of  neuron,  as  show  in  Figure  1.  This  3D 
geometry  can  be  further  meshed  in  the  finite  element 
model,  as  shown  in  Figure  2.  In  FEM  model,  the  neuron 
can  be  divided  into  three  different  regions:  main  cell  body, 
neurite  and  axon,  which  have  the  different  material 
properties.  In  the  current  work,  all  components  of  the 
neuron  are  modeled  as  both  the  elastic  and  the 


viscoelastic  materials.  The  standard  linear  model  (SLS)  is 
selected  for  the  viscoelastic  materials.  Based  on  the 
scanning  probe  indentation,  the  mechanical  properties  of 
the  neuron  was  reported  from  480  to  1600Pa  in  the 
frequency  range  of  30~200Hz.[3]  These  values  are 
considered  to  be  the  properties  of  the  cell  main  body  of  a 
neuron  since  the  indentation  test  is  usually  applied  on  the 
cell  main  body.  The  mechanical  property  ratio  of  the  axon 
and  the  neurite  to  the  cell  main  body  can  be  firstly 
assumed  to  be  one  and  then  gradually  increase  or  decrease 
to  match  their  simulated  deformation  field  with  the 
experimental  counterpart  (confocal  images  of  the 
deformed  neuron). 
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FIGURE  1. 2D  confocal  image  stack  of  an  isolated  neuron 
(0.5pm  between  each  image  from  the  top  surface  to  the 
substrate) . 


FIGURE  2.  (a)  3D  geometry  of  an  isolated  neuron  created 
from  2D  confocal  image  stack;  (b)  FEM  mesh  of  3D 
neuron  structure. 


simulated  deformation  field  and  the  experimental 
counterpart. 

CONCLUSIONS 

The  numerical  simulation  based  on  FEM  can  be  used  to 
determine  the  relative  strength  of  the  components  of 
neurons:  cell  main  body,  axon  and  neurite  by  matching 
the  simulated  deformation  field  of  a  neuron  and  its 
experimental  counterpart.  Then,  with  those  determined 
mechanical  properties,  the  neuron  local  mechanical 
response  and  the  relationship  between  the  local  response 
and  the  global  applied  deformation  can  be  obtained, 
which  will  be  the  first  step  of  building  the  mechanical 
injury  model  of  neuron. 
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EXPECTED  RESULTS 

In  FEM  simulations,  the  boundary  condition  of  the  neuron 
is  selected  similar  to  that  used  in  the  experiment.  In  the 
experiment,  the  neuron  is  assumed  to  be  perfectly  boned 
with  the  substrate,  and,  thus,  the  bottom  of  the  neuron  is 
clamped  in  the  simulation.  The  uniaxial  and  equibiaxial 
stretch  is  applied  on  the  neuron,  respectively.  The 
material  properties  of  the  neuron  components  are  the  set 
of  properties  which  give  the  least  difference  between  the 
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EXTENDED  ABSTRACT 


INTRODUCTION 

The  mechanism  of  mild  traumatic  brain  injury  (mTBI)  is 
directly  related  to  the  relationship  between  the  mechanical 
response  of  neurons  and  their  biological/chemical 
functions  since  the  neuron  is  the  main  functional 
component  of  brain.1  The  hypotheses  is  that  the  external 
mechanical  load  will  firstly  cause  the  mechanical 
deformation  of  neurons,  and  then,  when  the  mechanical 
deformation  of  neurons  reaches  to  a  critical  point  (the 
mechanical  deformation  threshold),  it  will  initiate  the 
chemical/biological  response  (e.g.  neuronal  function  loss). 
Therefore,  defining  and  measuring  the  mechanical 
deformation  threshold  for  the  neuronal  cell  injury  is  an 
important  first  step  to  understand  the  mechanism  of  mTBI. 

Typically,  the  mechanical  response  of  neurons  is 
investigated  based  on  the  deformation  of  in  vitro  model, 
in  which  the  neurons  are  cultured  on  the  elastic  substrate 
(e.g.  PDMS  membranes).  The  elastic  membrane  is 
deformed  by  the  external  load,  e.g.  equibiaxial  stretching. 
The  substrate  deformation  is  considered  to  be  the 
deformation  of  neurons  since  the  substrate  is  several 
orders  stiffer  than  the  neurons  and  the  neurons  are 
perfectly  bonded  with  the  substrate.  The  fluoresce 
method  is  typically  used  to  test  the  cell  injury,  e.g.  the  cell 
vitality  and  the  neuron  internal  ROS  level.1, 2 


EXPERIMENTAL  METHODS 

Neuroblastoma  cells  from  the  SH-SY5Y  cell  line  were 
used  to  measure  the  neuron  injury  threshold  caused  by  the 
mechanical  injury.  The  neurons  were  cultured  as 
monolayers  on  the  collagen-coated  six-well  FlexPlates 
(Flex,  Hillsborough,  NC,  USA)  in  DMEM  medium 
supplemented  with  10%  foetal  calf  serum  (FCS). 
Cultures  were  maintained  in  humidified  incubator  (5% 
C02,  37°C).  Cells  were  seeded  at  the  density  of  5000 
cells/cm2.  In  order  to  make  neurons  fully  differentiated, 
neurons  were  treated  with  10|iM  retinoic  acid  (RA)  over  4 
days  after  24hs  cell  culturing. 

The  cells  were  injured  using  Flexcell  5000  tension  system 
(Flex,  Hillsborough,  NC,  USA).  The  bottom  of  the 
FlexPlate  wells  are  made  by  PDMS  membrane  with  the 
thickness  of  0.2mm.  The  membrane  is  stretched  by  a 
negative  pressure  pulse.  The  schematic  of  the  cell 
stretching  is  shown  in  Figure  1.  The  cell  deformation  is 
same  as  the  deformation  of  the  PDMS  membrane.  In  the 
neuron  injury  evaluation  tests,  initially  an  arbitrarily  given 
deformation  e.g.  30%  will  be  applied  on  the  cultured 
neuron  monolayer,  which  should  be  strong  enough  to  kill 
most  of  the  PPDMS  membrane.  Then,  the  magnitude  of 
applied  strain  will  gradually  be  decreased  until  the  injury 
thresholds  are  measured.  The  loading  rate  is  about  0.01/s. 

The  nonlethal  injury  of  neurons  is  interested  since  it 
might  be  directly  related  to  the  mTBI.  The  nonlethal 


injury  can  be  defined  as  follows.  The  mild  injury  means 
that  the  neuron  can  still  survive  for  a  period  of  time  (24h) 
after  unloading.  The  critical  injury  strain  (£cr)  is  defined  as 
the  value  of  the  applied  external  deformation  beyond 
which  the  neuron  will  be  mildly  injured,  while  below 
which  there  will  be  no  detectable  injury.  Since  the  neuron 
injury  is  determined  based  on  the  statistical  results  in 
current  work,  the  applied  deformation  will  cause  about 
80%  of  the  neurons  in  each  test  (more  than  200  neurons 
are  counted)  to  have  the  detectable  injury  (nonlethal) 
defined  to  be  the  critical  injury  strain  (£cr).  The  upper 
limit  of  the  nonlethal  injury  of  neurons  (£^ax)  is  also 
defined.  If  the  applied  external  deformation  is  beyond  £^ax, 
the  neuron  will  finally  die.  £^ax  is  defined  as  the  value  of 
the  applied  external  deformation,  causing  about  30%  of 
neurons  to  die  in  each  test  (more  than  200  neurons  are 
counted). 

Dual  fluorescence,  live  cell  time-lapse  confocal 
microscopy  is  carried  out  for  identifying  the  neuron  injury. 
Two  different  fluorescent  reagents  will  be  used:  1)  a 
fluorescent  ROS  indicator,  2,7-dichlorofluorescin  (DCF,  1 
uM);  2)  a  HEPES-MEM  containing  0.5  uM  propidium 
iodide  (PI,  viability  indicator);  The  stretched  neurons  will 
be  loaded  with  both  regents.  An  increase  in  DCF  green 
fluorescence  strength  (usually  in  the  cell  cytoplasm) 
indicates  the  ROS  reactions  occurring  as  a  result  of  a 
higher  stress  level  of  the  cells  in  response  to  the  external 
load  (mechanical  stretch  and  contraction).  When  cell 
viability  changes  as  a  result  of  alteration  of  membrane 
permeability  in  response  to  the  mechanical  stretch  or 
contraction,  the  PI  will  enter  the  cells  and  binding 
DNAs/RNAs  (proteins  inside  nucleus),  resulting  in  red 
fluorescent  nucleus  structure. 


FIGURE  1.  Schematic  of  neuron  equibiaxial  stretching. 

RESULTS 

Figure  2  shows  the  florescence  images  of  SHSY5 
neuronal  cells.  The  cell  density  variation  caused  by  the 
stretching  is  measured  using  the  software  Image  J  (NIH). 
Figure  3  shows  the  neuron  density  variation  with  e  =10%, 
20%  strain  after  24h,  and  the  control  unit  after  24h  is  also 
shown  in  Figure  3.  The  results  show  that  the  equibiaxial 
strain  of  10%  is  the  upper  limit  of  the  non-lethal  injury 
(^maxX  which  caused  the  death  of  about  30%  neurons.  The 
lower  limit  of  the  non-lethal  injury  (£cr)  is  about  5%, 
which  is  measured  using  the  strength  of  DCF.  The  results 
of  each  well  are  based  on  the  average  values  of  15  images. 


FIGURE  3.  Fluorescence  images  of  neuron  density 
change,  (a)  control  unit;  (b)  £=10%;  (c)  £=20%. 


CONCLUSIONS 

The  nonlethal  neuron  injury  is  defined  and  measured 
based  on  the  In-Viro  cell  model  under  the  equibiaxial 
stretching.  For  SHSY5  neuronal  cells,  the  upper  limit  of 
the  non-lethal  injury  (£^ax),  which  caused  the  death  of 
about  30%  neurons.  The  lower  limit  of  the  non-lethal 
injury  (£cr)  is  about  5%, 
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Viscoelastic  mechanical  properties  of  biological  cells  are  commonly  measured  using  atomic  force  micro¬ 
scope  (AFM)  dynamic  indentation  with  spherical  tips.  A  semiempirical  analysis  based  on  numerical  simulation 
is  built  to  determine  the  cell  mechanical  properties.  It  is  shown  that  the  existing  analysis  cannot  reflect  the 
accurate  values  of  cell  elastic/dynamic  modulus  due  to  the  effects  of  substrate,  indenter  tip  size,  and  cell  size. 
Among  these  factors,  substrate  not  only  increases  the  true  contact  radius  but  also  interferes  the  indentation 
stress  field,  which  can  cause  the  overestimation  of  cell  moduli.  Typically,  the  substrate  effect  is  much  stronger 
than  the  other  two  influences  in  cell  indentation;  and,  thus,  the  cell  modulii  are  usually  overestimated.  It  is 
estimated  that  the  moduli  can  be  overestimated  by  as  high  as  over  200%  using  the  existing  analysis.  In  order 
to  obtain  the  accurate  properties  of  cells,  correction  factors  that  account  for  these  effects  are  required  in  the 
existing  analysis. 

DOI:  10.11 03/Phy sRe vE. 8 1 . 02 1 924  PACS  number(s):  87.10.Kn 


I.  INTRODUCTION 

There  is  a  close  relationship  between  the  physical  func¬ 
tions  of  cells  and  their  mechanical  behavior.  For  example, 
Cross  et  al.  recently  reported  that  the  stiffness  of  metastatic 
cancer  cells  is  70%  lower  than  that  of  healthy  cells  [1];  the 
heart  muscle  cells  loss  their  contractility  will  cause  the  heart 
failure  [2];  the  traumatic  brain  injury  (TBI)  might  be  caused 
by  the  neuron  death  under  the  severe  stretch  [3].  Cells  are 
typically  considered  as  continuum  materials  and  their  me¬ 
chanical  behavior  are  described  using  continuum  mechanics 
models.  The  parameters  associated  with  models  are  consid¬ 
ered  to  be  the  mechanical  properties  of  cells  which  can  be 
determined  using  experimental  measurements.  The  most 
commonly  used  models  are  elastic  and  viscoelastic  models 
based  on  which  the  elastic  or  dynamic  modulus  can  be  de¬ 
termined  [4-15]. 

AFM  is  considered  to  be  one  of  the  best  candidates  for 
accurately  measuring  the  load-displacement  relationship  of 
supersoft  materials  (~kPa)  at  a  small  scale.  Thus,  AFM  is 
becoming  one  of  the  most  popular  ways  to  measure  the  cell 
properties  [5-18].  Although  the  commercial  nanoindenter 
has  been  widely  used  to  measure  the  mechanical  properties 
of  materials  at  nanoscale  or  microscale,  the  applied  force 
range  of  nanoindenter  is  still  too  high  to  accurately  measure 
cell  properties.  For  example,  the  applied  force  range  for  the 
Hysitron  nanoindenter  is  from  30  nN  to  10  N  according  to 
the  menu  [19]. 

The  standard  sharp  AFM  tip  is  likely  to  cause  a  very  high 
stress  concentration  leading  to  highly  nonlinear  behavior. 
Also,  the  sharp  tip  may  even  tear  the  cell  membrane.  To  avert 
these  problems,  a  spherical  tip  is  usually  used  instead  of 
conventional  tip  in  cell  indentation  [7,18].  In  experiments, 
the  radius  of  the  indenter  tip  is  usually  chosen  under  40  jam 
[5-7,9,18],  compared  to  cells  which  are  typically  below 
40  jam. 
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The  cell  is  typically  mounted  on  a  substrate.  Since  the 
substrate  is  orders  of  magnitude  stiffer  than  the  cell  by  itself, 
the  indentation  force-displacement  (P-  S)  relationship  will  be 
influenced,  especially  for  low  cell  thickness  with  deep  inden¬ 
tation  [5,7].  Though  the  substrate  effect  is  minimal  in  shal¬ 
low  indentations,  it  is  not  practical  to  use  shallow  indentation 
to  cells  for  the  following  reasons.  In  shallow  indentations,  it 
is  highly  difficult  to  accurately  measure  the  contact  area  due 
to  the  surface  roughness  and  very  low  stiffness  of  cells  as 
well  as  the  presence  of  adhesive  force  between  the  indenter 
tip  and  cell  surface  [10,20,21].  It  is  very  important  to  under¬ 
stand  and  remove  the  substrate  effect  from  the  results  of  the 
AFM  indentation  measurement  in  order  to  obtain  the  true 
cell  properties. 

The  substrate  stiffening  effect  on  the  thin  film  indentation 
behavior  has  been  widely  investigated  [22-24].  Based  on  the 
approaches  initiated  by  Chen  [23]  and  Tu  et  al.  [24],  Dimi- 
triadis  et  al.  introduce  the  correction  terms  into  the  Hertz 
contact  model  to  correct  the  substrate  stiffening  effect  on  the 
elastic  modulus  of  cell  [25].  Mahaffy  et  al.  [7]  developed  an 
analysis  to  obtain  the  viscoelastic  properties  of  cells  attached 
on  the  rigid  substrate  from  the  AFM  dynamic  indentation 
measurement.  In  this  approach,  the  solution  to  the  problem 
of  cell  indentation  with  substrate  was  assumed  to  be  a  series 
expansion  of  the  standard  Hertz  elastic  contact  solution  ex¬ 
tended  to  linear  viscoelastic  materials.  All  of  the  above  re¬ 
sults  and  other  reported  cell  properties  measured  using  AFM 
indentation  [6,9]  are  still  based  on  the  Hertz  contact  solution 
which  is  the  solution  derived  for  the  semi-infinite  elastic  con¬ 
tact  problem.  In  addition,  the  geometric  characteristics  of  the 
cell  (the  cell  diameter  and  thickness  as  well  as  the  indenter 
tip  size)  are  not  considered,  which  may  significantly  affect 
the  cell  indentation  response.  For  example,  the  substrate 
stiffening  effect  may  not  be  just  a  function  of  cell  thickness 
but  also  dependent  on  the  indenter  tip  size  and  cell  diameter. 
Therefore,  it  is  highly  necessary  to  establish  an  effective 
analysis  to  identify  these  effects  and  show  the  intrinsic  prop¬ 
erties  of  cell. 

In  this  paper,  the  cell  mechanical  behavior  is  investigated 
using  the  semiempirical  approach  based  on  numerical  simu¬ 
lations.  In  the  simulations,  the  relationship  between  contact 
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FIG.  1.  (Color  online)  The  schematic  of  cell  indentation  and 
simplified  disklike  shape. 

radius  and  indentation  displacement  is  established  and  the 
substrate  effect  is  identified.  The  effects  of  indenter  tip  size, 
cell  thickness,  and  cell  radius  on  the  indentation  behavior  are 
examined.  The  correcting  factors  are  introduced  into  the  ex¬ 
isting  indentation  analysis;  then  cell  properties  can  be  more 
accurately  determined.  This  study  can  help  us  to  understand 
the  intrinsic  mechanical  properties  of  cells  and  can  provide  a 
useful  guideline  for  building  the  relationship  between  me¬ 
chanical  behaviors  and  biological  functions. 


is  represented  by  25  000  four-node  axisymmetric  elements 
with  reduced  integration.  The  indenter  tip  radius  (. R )  is  se¬ 
lected  from  1  —  15  n m,  which  is  commonly  used  in  cell  in¬ 
dentation  [5,7].  The  different  cell  radii  (r)  are  selected  as  r 
=  10,  20  yarn  (typically  less  than  20  yam  in  radius).  The  dif¬ 
ferent  thicknesses  (h)  are  selected  as:  h  =  5,  10  yam  (typically 
larger  than  3  yam).  In  order  to  obtain  an  accurate  contact 
radius,  the  size  of  surface  elements  in  the  contact  area  is  set 
to  less  than  0.5%  of  the  indenter  tip  radius.  The  Prony  series 
coefficients  of  SLS  model  are  g  =  0.3  —  0.9  and  r=l  — 10  s 
as  well  as  the  instantaneous  modulus  range  is  Einst3Lnt=10 
— 100  kPa  in  FEM  simulations.  A11  degrees  of  freedom  of 
nodes  on  the  bottom  of  the  cell  are  constrained  to  simulate 
the  condition  that  the  cell  is  fully  adhered  to  the  rigid  sub¬ 
strate  surface.  A11  indentation  simulations  are  performed 
based  on  displacement  control.  In  order  to  examine  the  true 
substrate  effect,  the  indentation  displacement  selected  in  the 
simulations  coincide  with  the  lowest  experimentally  applied 
value  [5,7]. 


II.  COMPUTATIONAL  METHODS 

In  the  present  work,  the  cell  indentation  behavior  is  stud¬ 
ied  using  numerical  simulations  based  on  finite  element 
modeling  (FEM).  After  cells  mounted  on  substrate,  cells  will 
expand  themselves  and  the  size  along  the  normal  direction  of 
substrate  is  much  less  than  the  lateral  size  (as  shown  in  Fig. 
1).  Cell  geometry  is  typically  described  using  two  param¬ 
eters:  cell  average  radius  r  (along  the  lateral  direction)  and 
thickness  h  (along  the  normal  direction  of  substrate) 
[4,6,10,11,26-29]. 

The  cell  is  usually  modeled  as  a  continuous  isotropic  lin¬ 
ear  viscoelastic  material.  The  standard  linear  solid  (SLS) 
model  is  used  to  determine  the  dynamic  modulus  of  cells  as 
shown  in  Fig.  2.  In  the  SLS  model,  the  relaxation  modulus  is 
given  by 


E{t)  =  Ei  +  E2e~tlT,  (1) 

where  Ex  and  E2  are  elastic  modulus  of  springs  and  r  is  the 
relaxation  time,  r=  rj/E 2,  where  77  is  the  viscosity.  The  sub¬ 
strate  and  indenter  tip  are  assumed  to  be  rigid  since  they  are 
usually  several  orders  stiffer  than  the  cell. 

A11  FEM  simulations  are  performed  using  commercial 
code  ABAQUS  v.6.8.  Cells  are  geometrically  simplified  as 
two-dimensional  (2D)  axisymmetric  disks.  The  spherical  in¬ 
denter  tip  is  modeled  as  a  2D  axisymmetric  surface.  The  cell 
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FIG.  2.  (Color  online)  The  standard  linear  solid  model  of  vis¬ 
coelastic  material. 


III.  CONTACT  RADIUS  IN  CELL  INDENTATION 

In  indentation  tests,  the  accuracy  of  the  result  is  strongly 
influenced  by  the  accuracy  demonstrated  in  obtaining  the 
indentation  contact  radius.  In  the  following  subsections,  we 
will  discuss  the  effects  of  geometric  parameters  (indenter  tip 
radius  R ,  cell  thickness  h  and  cell  radius  r)  on  the  indentation 
contact  radius  in  cell  indentation  and  the  corresponding 
physical  mechanisms. 

Based  on  dimensional  analysis,  the  normalized  contact 
radius  can  be  described  as  a  function  of  the  normalized  in¬ 
dentation  displacement  (S/h),  indenter  tip  radius  (Rlh),  cell 
radius  (r/R),  and  indentation  force  (PI Eh2): 

a/R  =  F(  S/h ,  R/h,  r/R,  PI  Eh2) . 

The  related  geometric  parameters  are  R,  h,  and  r.  Although  it 
is  reported  recently  that  cells  might  be  compressible  [30,31], 
cells  are  commonly  assumed  as  incompressible  (^—0.5).  In 
the  present  work,  the  incompressible  assumption  is  still  used, 
and  thus,  v  is  not  considered  as  a  variable  in  the  present 
paper. 


A.  Effects  of  geometric  parameters  on  contact  radius 

The  variation  in  anum/R  with  S/h  is  shown  in  Figs. 
3(a)-3(d),  where  anum  is  the  numerical  solution  of  the  contact 
radius.  In  the  figure,  the  indenter  tip  radius  R=  1,  15  yam,  the 
cell  radius  r=  10,  20  yam  and  the  cell  thickness  h  =  5,  10  yam 
respectively.  The  FEM  results  are  shown  as  thin  lines  (not 
smooth).  For  the  sake  of  reference,  the  normalized  Hertz 
contact  radius,  ah/R,  is  shown  as  the  dashed  line.  It  is  seen 
that  when  both  R  and  h  are  small,  the  numerical  solution 
( anum )  is  very  close  to  the  Hertz  solution  (ah)  and  anum  is  not 
sensitive  to  r  [the  FEM  result  of  the  case  with  r=  10  yam  is 
overlapped  with  the  Hertz  solution  in  Fig.  3(a)].  With  the 
increase  in  R  or  h,  anum  deviates  from  ah  and  the  deviation 
increases  with  S/h.  When  both  R  and  h  are  large,  anum  in¬ 
creases  with  r.  The  effects  of  R,  h,  and  r  on  anum  are  coupled 
with  each  other.  In  addition,  for  a  smaller  R  and  a  larger  h. 
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FIG.  3.  (Color  online)  The  relationship  between  normalized 
contact  radius,  a/R ,  and  normalized  indentation  displacement,  81  h. 

anum<ahi  f°r  a  larger  R  and  a  smaller  h ,  anum>ah. 

For  all  cases  in  the  present  study  (1<R<15  fi m,  h 
=  5,10  /z m,  r=  10,20  /zm,  S/h<  10%),  the  numerical  solu¬ 
tion  of  contact  radius  can  be  fitted  as 


Normalized  IndenterTip  Radius,  R/h 
FIG.  4.  (Color  online)  The  fitting  parameter  k  in  Eq.  (2). 
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(2) 


where  k=k(R,h )  is  the  fitting  parameter,  which  also  depends 
on  the  cell  radius,  r.  The  fitting  curves  of  FEM  results  are 
shown  as  thick  solid  lines  (smooth)  in  Figs.  3(a)-3(d).  For  all 
cases  except  for  R=l  /zm  and  8=1  /z m,  k  can  be  fitted  as 
k=ki  +  k2R/h ,  as  shown  in  Fig.  4.  The  fitting  parameter,  kx,  is 
essentially  a  constant  and  kx~ -0.083.  The  fitting  parameter 
k2~0.48  for  r= 20  /zm.  However,  k2  decreases  with  the  in¬ 
crease  in  8  or  h  when  r  reduces  to  10  /zm,  as  shown  in  the 
Fig.  4. 


B.  Basic  mechanisms 

It  is  evident  that  all  geometric  parameters  R,  h,  and  r  will 
influence  the  contact  radius  due  to  nonlinear  geometry,  sub¬ 
strate,  and  boundary  effects.  Substitution  of  the  expression  of 
k-kx  +k2R/h  into  Eq.  (2)  leads  to 

anum  —  (RS)112  +  kx8+  k2—S.  (3) 

h 

If  8<R  and  8<h,  Eq.  (3)  converges  to  the  Hertz  solution. 
The  second  term,  kxS,  in  Eq.  (3)  is  the  geometric  nonlinear 
correction  term  of  the  contact  radius  when  the  condition  8 
<  R  is  not  satisfied.  The  geometric  nonlinear  effect  increases 
with  the  increase  in  8  or  the  decrease  in  R.  In  the  present 
paper,  the  inclusion  of  the  first-order  term  of  8  is  accurate 
enough  to  describe  the  geometric  nonlinear  effect.  In  fact,  for 
very  small  R  and  large  8 ,  higher-order  terms  of  8  are  re¬ 
quired.  For  example,  at  R=  l  /zm  and  £=1  /zm  the  first- 
order  term  of  8  is  not  accurate  enough  (Fig.  4).  Since  kx  is 
negative,  the  nonlinear  geometric  effect  will  reduce  the  con¬ 
tact  radius  compared  with  the  Hertz  solution.  Generally,  the 
nonlinear  geometric  effect  can  be  reduced  with  a  larger  in- 
denter  tip  size. 
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FIG.  5.  (Color  online)  The  schematics  of  cell  indentation  with 
and  without  substrate. 

The  third  term  in  Eq.  (3)  can  be  considered  as  the  sub¬ 
strate  effect.  Figure  5  shows  the  schematics  of  cell  indenta¬ 
tion  with  and  without  substrate.  The  cell  compression  in  the 
thickness  direction  (introduced  by  the  indentation)  produces 
the  cell  stretching  in  the  radial  direction  due  to  Poisson’s 
effect.  Both  the  compression  and  the  stretching  cause  the 
penetration  depth,  8in,  to  be  less  than  the  displacement  of 
indenter  tip,  8.  For  example,  £ ~  2  Sin  in  the  Hertz  solution.  If 
a  cell  is  fully  adhered  to  a  substrate,  the  substrate  will  con¬ 
strain  the  bottom  surface  of  the  cell.  Due  to  the  small  cell 
thickness,  this  constraint  will  reduce  the  cell  stretching  in  the 
radial  direction.  This  constraint  effect  increases  the  penetra¬ 
tion  depth,  8in,  under  the  same  applied  8  compared  to  the 
case  without  the  substrate.  Thus,  the  contact  radius  will  be 
larger  than  the  Hertz  contact  radius  due  to  the  substrate  ef¬ 
fect.  The  substrate  effect  on  the  contact  radius  increases  with 
the  decrease  in  h  or  increases  with  R. 

The  boundary  effect  reduces  the  contact  radius  by  affect¬ 
ing  the  extent  of  the  cell  deformation  along  the  radial  direc¬ 
tion  caused  by  Poisson’s  effect.  For  given  R ,  h ,  and  8 ,  a 
smaller  cell  will  deform  easier  than  a  larger  cell  along  the 
radial  direction  due  to  the  interaction  between  the  free  sur¬ 
face  and  the  indentation  stress  field.  This  can  be  further  sup¬ 
ported  by  comparing  the  displacements  along  the  radial  di¬ 
rection,  un ,  of  both  cells.  The  results  show  that  un  on  the 
free  surface  of  the  smaller  cell  is  much  larger  than  un  at  the 
same  radial  position  of  the  larger  cell.  The  easier  lateral  de¬ 
formation  will  reduce  the  penetration  depth  and  thus  reduce 
the  contact  radius.  This  effect  increases  with  8  or  R  or  de¬ 
creases  with  the  increase  in  r. 

In  summary,  compared  with  the  Hertz  contact  solution, 
the  nonlinear  geometry  and  the  boundary  effects  reduce  the 
contact  radius,  while  the  substrate  effect  increases  the  con¬ 
tact  radius.  With  a  large  R/h  (R/h> 0.2),  the  substrate  effect 
is  stronger  than  the  nonlinear  geometric  and  the  boundary 
effects;  and,  thus,  the  contact  radius  is  underestimated  using 
Hertz  contact  radius.  With  a  small  R/h  (R/h< 0.2),  the  sub¬ 
strate  effect  is  smaller  than  the  nonlinear  geometrical  effect 
(the  boundary  effect  is  very  weak  under  this  condition  since 
R/r  is  very  small);  and,  thus,  the  contact  radius  is  overesti¬ 
mated  using  Hertz  contact  radius.  Since  typically  R/h> 0.2 
in  the  cell  indentation,  the  contact  radius  is  underestimated 
using  the  Hertz  contact  radius.  This  underestimation  is 
smaller  for  a  smaller  R/h.  In  addition,  the  relationship  be¬ 
tween  contact  radius  and  indentation  displacement  is  only 


FIG.  6.  (Color  online)  The  relationship  between  the  normalized 
indentation  force  and  the  normalized  indentation  displacement. 


dependent  on  the  geometric  parameters  but  not  the  material 
properties  of  cell  (. E ,  g,  and  r).  This  has  been  validated  using 
FEM  simulations  based  on  the  different  values  of  E,  g,  and  r 
(£=10-100  kPa,  g= 0.3  — 0.9,  and  r=l-10  s). 


IV.  EFFECT  OF  CONTACT  RADIUS  ON  ELASTIC 
MODULUS 

Since  the  Hertz  solution  cannot  accurately  describe  the 
contact  radius  in  cell  indentation,  in  this  section,  we  will 
discuss  the  effect  of  the  Hertz  contact  radius  on  determining 
the  cell  modulus  based  on  both  elastic  and  viscoelastic  ma¬ 
terial  models.  In  order  to  obtain  the  correct  cell  modulus,  the 
correcting  factor  for  the  Hertz  contact  radius  is  required. 


A.  Elastic  materials 


For  a  rigid  spherical  indenter,  the  indentation  force  is 
given  by  [32] 


4  E  a3 

P= - , 

31  -v1  R 


(4) 


where  v  is  the  Poisson’s  ratio  and  E  is  the  elastic  modulus.  If 
the  deformation  is  infinitesimal,  8<R ,  the  contact  radius  is 
commonly  approximated  as  the  Hertz  solution:  a~ah=  ^ Jrs . 
Thus,  the  indentation  force  can  be  simplified  as  [32] 

4  EiRSm 


The  cell  elastic  modulus  can  be  directly  calculated  from  the 
measured  Ph~  8 relationship  based  on  Eq.  (5).  One  example 
of  the  indentation  P  —  8  curve  computed  from  the  numerical 
simulation  is  shown  in  Fig.  6.  Pnum  is  the  numerical  solution 
of  indentation  force.  The  indentation  force  Pc  calculated 
from  Eq.  (4)  based  on  anum  and  the  Hertz  solution  of  inden¬ 
tation  force  Ph  calculated  from  Eq.  (5)  are  also  indicated  as  a 
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FIG.  7.  (Color  online)  The  correction  factors  to  compensate  for 
the  effect  from  the  Hertz  contact  radius  on  the  cell  elastic  modulus, 

/cl  and  fc2. 


dashed  line  and  a  dash-dot  line,  respectively.  It  can  be  seen 
that  Pnum>Pc>Ph •  Pnum>Pfi  even  at  a  small  8  and  Pc 
~  Pnum  when  S/h<  4%.  With  the  increase  in  8 ,  Pnum  begins 
to  deviate  from  Pc.  This  deviation  is  caused  by  the  stress 
stiffening  effect  from  substrate,  which  will  be  discussed  in 
detail  in  the  next  section.  The  difference  between  Pc  and  Ph 
results  from  the  inaccurate  contact  radius  by  the  Hertz  solu¬ 
tion.  Based  on  Eq.  (4),  the  underestimation  of  contact  radius 
will  cause  the  overestimation  of  elastic  modulus  by  a  factor 
fci  =  (anUJah, )3  as  shown  in  Fig.  7. 

The  more  popular  way  to  determine  the  elastic  modulus  in 
nanoindentation  tests  is  based  on  the  contact  stiffness  dP Id 8, 
which  can  be  calculated  by  the  first  derivative  of  P  with 
respect  to  8  in  Eq.  (5).  The  elastic  modulus  can  be  expressed 
as  [33] 

E  dP  1 
(1  -  i?)  d82ah  ’ 


where  the  contact  stiffness  dP/d8  can  be  measured  from  the 
P-8  curve  in  experiments.  Since  dP  I  d8  is  easily  measured  in 
experiments  (from  the  initial  stage  of  the  unloading  curve 
slop),  this  way  is  more  popular  to  measure  the  cell  elastic 
modulus.  However,  if  the  contact  radius  is  described  by  Eq. 
(2),  Eq.  (5)  will  change  to 


_  4 _E_(jRS+kS)^ 
P~  3  1  -  v2  R 


(7) 


where  n i  =  1 ,  n2=3,  n3  =  3,  n4=  1.  When  m- 1,  i.e.,  only 
the  first  term  of  the  sequence  in  Eq.  (8)  is  considered,  then  it 
reduces  to  Eq.  (6).  Thus,  if  anum>ah ,  then  Eq.  (6)  can  over¬ 
estimate  the  elastic  modulus  by  a  factor  fc2  in  quasistatic 
indentation: 


fc2  ~ 


2 m 

3  ah  ’ 


(9) 


4  +2 

£(<?)  =  2  £m.  L  =  nm  — —  (k)m~l(R)l~m,2(S)ml2 . 

m=  1  ^ 

(10) 

fC2  is  also  shown  in  Fig.  7.  It  can  be  seen  that /cl  >fc 2  for  a 
given  8.  Thus,  Although  Eqs.  (5)  and  (6)  are  equivalent  when 
a  ~ ah=^R8,  Eq.  (6)  will  cause  a  large  overestimation  for  the 
elastic  modulus  than  Eq.  (5)  when  the  contact  radius  is  de¬ 
scribed  using  Eq.  (2).  In  addition,  E  is  not  a  constant  but  a 
function  of  the  indentation  displacement  8. 


B.  Viscoelastic  materials 

During  the  loading  stage,  the  indentation  force  of  vis¬ 
coelastic  material  can  be  calculated  using  the  correspondence 
principle:  replacing  the  time-independent  constant  in  Eq.  (5) 
by  the  corresponding  differential  operators  of  the  viscoelastic 
constitutive  model  [34]. 


p(t)= 


3(1  -  i?)R 


f 

Jo 


ds 


(id 


In  dynamic  indentation  with  displacement  control,  the  ap¬ 
plied  displacement  profile  is  8(t)  =  8q+ A£sin(cor).  The  in¬ 
dentation  force  can  be  also  considered  as  a  direct  indentation 
force  component  superimposed  with  an  oscillatory  compo¬ 
nent  P(t)  =  P0+AP  sin (tot+ffi).  The  storage  modulus  and  the 
loss  modulus  can  be  commonly  determined  based  on  the 
contact  stiffness  (A/VA£)  (Ref.  [35]): 


E'  1  A  P 

- 7  = - — :  cos  cb, 

1  -  v1  2 ah  A  8 


(12) 


E”  1  A  P 

- «  = - - — sin  (/). 

1  -  i?  2 ah  A  8 


(13) 


However,  when  the  contact  radius  is  described  by  Eq.  (2), 
based  on  Eq.  (7),  Eq.  (11)  will  change  to 


P(f)  = 


3(1  -  i/)R 


f 

Jo 


d{^R8(t  -  s)  +  k8(t  -  s)}3 

E(s) - ds . 

o  ds 


(14) 


Then, 


Based  on  the  first  derivative  of  P  with  respect  to  8  in  Eq.  (7), 
the  elastic  modulus  is  given  by 
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Pit)  =  2  nm{k)m-\R)x~ml2  f  E(s )—  - ds, 

3 1  -  i/m=l  J0  ds 


(15) 


where 


ml  2+1 


§(t)ma+l  =  (<50  +  A<5  sin  cot) 

=  ^/2+1  +  <5jf2A S sin  cot  +  o(AS), 
m  =  1—4. 


(16) 


Let 


TVl  “I-  2 

&s0)  =  X  nm—(kr-1(R)1-m,2(80r'2.  (17) 

vn—  1  ^ 


Therefore,  the  oscillatory  part  is  given  by 


A P  sin(w?  +  cf>)  =  ^y*-^£(S0)[E'(co)2 

+  £’"(ft>)2]1/2A<5sin((wf  +  cf>),  (18) 


FIG.  8.  (Color  online)  The  relationship  of  the  overall  indentation 
correction  factor  ftotl  and  normalized  indentation  depth,  S/h. 

A.  Elastic  materials 


tan  cf)  =  E"(co)/Ef(co), 


E'(gj)  3  A P  1 
1-v2'4A^(^)COS  (f>’ 


E"(<o)  3  A  P  1 

- ~  = - ; — -sin  cb. 

l-v1  4  A  <5  £(<50)  V 


(19)  The  overall  influences  of  both  the  underestimation  of  con¬ 
tact  area  and  the  stress  stiffening  effect  from  the  substrate 
can  be  simply  represented  by  a  factor  ftot\=Pnum/Ph  based 

(20)  on  Eq.  (4),  which  are  shown  as  thin  lines  in  Fig.  8.  It  can  be 
seen  that  for  all  values  of  R ,  ftotX  ~  1  at  a  very  small  S/h;  and 
ftotX  increases  with  S/h  and  R.  It  shows  that  the  elastic  modu¬ 
lus  is  overestimated  by  a  factor,  ftotl,  based  on  Eq.  (4)  in  the 

^  '  quasistatic  indentation.  The  factor  includes  two  components: 


Similar  to  the  analysis  for  elastic  materials,  if  a  >ah 
=  ^RS,  the  conventional  analysis  [Eqs.  (12)  and  (13)]  will 
overestimate  the  complex  modulus  with  the  same  factor: 


2  g(4) 

3  a  hi  do) 


(22) 


From  Eqs.  (20)  and  (21),  the  dynamic  modulus  components 
are  not  only  functions  of  the  frequency  of  the  oscillatory 
load,  a ),  but  also  functions  of  direct  indentation  displace¬ 
ment,  S0.  The  phase  difference,  </>,  between  the  indentation 
force  and  the  indentation  displacement  is  not  affected  by  the 
correction  of  contact  radius.  The  above  results  have  been 
validated  by  FEM  simulations  with  the  SLS  model  (g  =  0.3 
—  0.9,  T—  1  10  s,  and  £’instant=  10—  100  kPa)  and  actually 

fc 2  does  not  depend  on  the  parameters  in  SLS  model. 


V.  NONLINEAR  STRESS  STIFFENING  EFFECT 


fwt\=fs+fc\- 1.  (23) 

where  fs  is  from  the  substrate  stiffening  effect  and  /cl  is  from 
the  underestimation  of  contact  radius.  For  example,  fcl 
=  1.92,  ftot\= 2.3  at  S/h=  10%,  and,  thus,  fs=  1.38  for  the  cell 
with  h-  5  yam,  R  =15  /z m,  and  r=20  jmm.  This  means  that 
the  overall  overestimation  of  the  elastic  modulus  is  130%, 
wherein  the  underestimation  of  the  contact  radius  contributes 
to  about  92%,  and  the  rest  from  the  substrate  stiffening  effect 
is  about  38%. 

Figure  9  shows  the  variation  in  ftotl  with  R/h  at  S/h 
=  10%.  ftotl  increases  with  R/h.  With  a  larger  cell  radius  r 
=  20  /zm,  ftotl  is  essentially  insensitive  to  h.  Further,  with  a 
smaller  thickness  h-  5  /zm,  ftotl  is  insensitive  to  r;  but  the 
effect  of  r  becomes  larger  with  h.  The  effect  of  r  on  ftotl 
arises  mainly  from  the  influence  of/cl,  which  can  be  reduced 
by  selecting  smaller  R/h. 

For  a  given  cell  radius,  r,  ftotl  can  be  fitted  as  a  function: 


In  this  section,  we  will  discuss  the  substrate  stiffening 
effect  on  the  cell  modulus  based  on  both  elastic  and  vis¬ 
coelastic  material  models.  This  effect  will  cause  a  significant 
overestimation  of  the  cell  elastic  or  dynamic  modulus  based 
on  the  Hertz  solution.  In  order  to  get  the  accurate  cell 
moduli,  the  correcting  factor  for  the  substrate  stiffening  ef¬ 
fect  is  required  for  the  conventional  analysis. 


ftotl  =  PnuJPh  =  1  +  b(S/h) 1/2  +  c(S/h) .  (24) 

The  fitting  results  of  the  cells  with  r=20  /zm  and  h 
=  5  /zm  are  shown  as  thick  lines  in  Fig.  8  (they  are  over¬ 
lapped  with  the  FEM  results).  The  coefficients  b  and  c  in  Eq. 
(24)  can  be  fitted  as  functions  of  R/ h\b=bl(R/ h)m+b2,  c 
=  cxR/ h  +  c2,  as  shown  in  Fig.  10.  For  cells  with  h  =  5  /zm, 
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(a)  The  Ratio  between  Indenter  Tip  Radius  and  Cell  Thickness,  R/h 


FIG.  9.  (Color  online)  The  comparison  between  the  overall  cor¬ 
rection  factors,  ftoti  and  ftot2  of  the  elastic  modulus. 


these  constants  can  be  found  as  Z^  —  0.73,  Z?2~0.27,  and  c1 
«3.0,  c2~-0.9. 

Pnum  can  be  expressed  as  that  Ph  multiplied  by  a  function: 
g(S/h)  =  1  +b(S/ h)ll2+c(S/ h),  which  represents  the  effect  of 
geometric  parameters  (small  R ,  large  S ,  and  low  h  with  sub¬ 
strate).  Thus, 


4  e/rS2'2  (  S\ 
~  3  1  -  v2  Ah)' 


(25) 


Then  based  on  the  contact  stiffness,  dPnum/dS ,  Eq.  (6)  will 
change  to 


P  3  dP 
_E_  =  2armrl 

1  -  v2  4  b 

3 

£(S)  =  /?1/2E  lj—^8n/2h(l-m)l2,  (27) 

m=  1  ^ 

where  /x  =  1 ,  /2  =  &,  /3  =  c.  Therefore,  if  a>ah=^R8  and 
with  the  substrate  effect,  The  more  popular  method  based  on 
the  contact  stiffness  [Eq.  (6)]  will  overestimate  the  elastic 
modulus  by  a  factor  ftot2: 

ftot2=^Jy  (28) 

3  ah{S) 

It  can  be  seen  that  ftot2>ftotb  as  shown  in  Fig.  9.  The  over¬ 
estimation  of  elastic  modulus  can  be  as  high  as  200%  based 
on  the  more  popular  way  [Eq.  (6)],  wherein  the  underestima¬ 
tion  of  a  contributes  about  130%,  as  shown  in  Fig.  7,  and  the 
rest  arising  from  the  stress  stiffening  effect. 


B.  Viscoelastic  materials 

For  the  viscoelastic  materials  with  substrate  effect,  the 
components  of  complex  modulus  are  given  by 


(b)  The  Ratio  between  Indenter  Tip  Radius  and  Cell  Thickness,  R/h 

FIG.  10.  (Color  online)  The  fitting  parameters  h  and  c  in  Eq. 
(24). 


E’(aj)  3  A P  1 
l-i/"4A^(4)C°S 

(29) 

E’{w)  3  A  P  1 

1-^'4A^(4)C°S 

(30) 

tan  </>  =  £"(«)/£'(«), 

(31) 

3 

where  ^(^0)=/?1/2  2  Similar  to  elastic  ma- 

m=  1 

terials,  if  a>ah=\lRS and  with  the  substrate  effect,  Eqs.  (12) 
and  (13)  will  overestimate  the  complex  modulus  by  the  same 
factor: 


021924-7 


GUOXIN  CAO  AND  NAMAS  CHANDRA 


PHYSICAL  REVIEW  E  81,  021924  (2010) 


_  2 gsp) 

3 ah(S0) ' 


(32) 


From  Eqs.  (29)-(31),  it  can  be  seen  that  the  phase  difference 
is  not  affected  by  the  substrate  effect;  and  the  components  of 
the  intrinsic  dynamic  modulus  can  be  provided  by  Eqs.  (29) 
and  (30).  The  above  results  have  been  validated  using  FEM 
simulations  with  the  SLS  model  (g  =  0.3~0.9,  r=l~10  s, 
and  Einstant  —10  ^  1 00  kPa)  and  actually  the  substrate  effect 
is  not  dependent  on  the  values  of  g  and  r  used  in  the  SLS 
model. 


VI.  SUMMARY  AND  CONCLUSIONS 

In  this  work,  a  semiempirical  analysis  is  built  based  on 
FEM  simulations  to  determine  the  cell  mechanical  properties 
from  the  spherical  AFM  indentation  response.  The  results 
show  that  the  existing  indentation  analysis  can  overestimate 
the  dynamic  modulus  by  as  high  as  200%  in  some  cases.  The 
overestimation  is  mainly  caused  by  the  inaccurate  determi¬ 
nation  of  the  indentation  contact  radius  and  the  improper 
accounting  of  the  substrate  effect.  The  effects  of  inaccurate 
contact  radius  and  substrate  depend  on  the  geometric  param¬ 
eters:  indentation  tip  radius,  cell  radius,  and  cell  thickness 
but  do  not  depend  on  the  material  properties  ( E ,  g ,  and  r). 
Correction  factors  have  been  proposed  in  this  work  that  re¬ 


late  the  geometric  parameters  to  the  mechanical  properties  in 
quasistatic  and  dynamic  indentation  tests. 

The  inaccurate  determination  of  the  contact  radius  is 
caused  by  the  nonlinear  geometric,  finite  boundary,  and  sub¬ 
strate  effects.  Substrate  not  only  affects  the  value  of  the  con¬ 
tact  radius  but  also  stiffens  the  cell.  The  substrate  effect  in¬ 
creases  with  an  increase  in  indenter  tip  radius,  indentation 
displacement,  or  a  decrease  in  cell  thickness.  The  nonlinear 
geometric  effect  increases  with  a  decrease  in  indenter  tip 
radius.  The  boundary  effect  reduces  with  cell  thickness  or/ 
and  indenter  tip  radius.  When  the  indenter  tip  radius  is  small, 
the  nonlinear  effect  is  dominant  compared  to  substrate  and 
boundary  effects.  When  the  tip  radius  is  large,  the  substrate 
effect  is  dominant.  To  remove  these  effects,  the  correction 
factors  ftoti  and  ftot2  are  induced,  composed  of  the  contact 
area  correction  factor  and  the  stress  stiffening  factor.  They 
increase  with  the  increase  inR/h  or  the  decrease  in  h.  For  a 
small  h ,  the  correction  factors  are  not  sensitive  to  r;  for  a 
large  h ,  they  decrease  with  r. 
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INTRODUCTION 

Blast  induced  traumatic  brain  injury  (bTBI)  is  signature 
injury  in  recent  combat  scenarios  involving  improvised 
explosive  devices  (IEDs).  The  exact  mechanisms  of  bTBI  are 
still  unclear  and  protective  role  of  helmet  and  body  armor  is 
often  questioned  [1-3].  High  Fidelity  finite  element  models 
involving  fluid  structure  interaction  are  built  in  order  to 
understand  effectiveness  of  helmet  in  mitigating  early  time 
blast  induced  mild  traumatic  brain  injury. 

METHODOLOGY 

Two  dimensional  plane  strain  finite  element  models  of 
helmet-head  under  shock  loading  are  studied  to  compare 
effectiveness  of  helmet.  Figure  1  shows  the  configuration  of 
setup. 


Figure  1  :  Simulation  setup.  Frontal  blast  is  simulated. 

Our  blast  scenarios  are  simulated  by  first  positioning  the 
head  model  in  an  atmosphere  of  air  at  ambient  conditions  as 
shown  in  Fig.l.  Shock  wave  is  generated  by  releasing  high 
pressure  compressed  air  into  atmospheric  air  at  time  equal  to 
zero.  The  pressure  and  thickness  of  compressed  air  domain  is 
selected  so  as  to  generate  nonlethal  blast  wave.  The  structure  of 
this  blast  wave  is  illustrated  in  Fig.  2. 
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Figure  2  :  Wave  form  of  approximated  air  blast  structure  of 
0.3  MPa  magnitude 

Our  head  model  consists  of  skull ,  facial  bones  ,  neck  bones 
and  brain.  We  have  used  same  material  properties  for  skull  , 
facial  bones  and  neck  bones  hence  we  do  not  distinguish  them 
as  separate  component,  skull  ,  facial  bones  and  neck  bones 
together  will  be  referred  as  ’skull’  henceforth.  The  geometry  of 
these  components  is  obtained  by  segmentation  of  MRI  dataset 
available  from  visible  human  project  of  the  National  Library  of 
Medicine  [4].  Since  we  are  interested  in  2-dimensioanl  head 
model  ,  the  central  slice  of  MRI  dataset  is  chosen.  These 
geometries  are  imported  into  finite  element  software  [5]  and 
then  meshed  to  generate  2-dimensional  plain  strain  finite 
element  head  model.  2-dimensioanl  geometries  were  selected 
so  that  the  analysis  would  not  be  overly  complex  and 
prohibitively  expensive. 

The  brain  tissue  is  modeled  as  linear,  isotropic,  viscoelastic 
material  with  properties  adopted  from  Taylor  et  al.[2].  Standard 
Linear  Solid  (SLS)  model  is  used  to  characterize  shear 
response.  The  skull  is  modeled  as  linear,  elastic,  isotropic 
materials  based  on  material  models  suggested  in  the  literature 
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[6-9].  The  Kevlar  helmet  is  modeled  as  hollow  hemiellipsoid 
with  a  constant  thickness  and  offset  from  the  skull  as  described 
by  Reynosa  [10],  with  transversely  isotropic  elastic  material 
properties  given  by  Aare  and  Keliven  [11].  Dry  air,  the  medium 
through  which  blast  wave  propagates  is  modeled  as  ideal  gas 
equation  of  state. 

RESULTS 

We  have  studied  how  helmets  influence  the  blast-induced 
mechanical  loads  in  the  brain.  Helmeted  and  non  helmeted 
response  is  compared  on  the  basis  of  pressure  in  the  helmet 
cranium  subspace  ,  contact  pressure  on  the  outer  surface  of  the 
skull  and  pressure  and  shear  stresses  (mises)  in  the  various 
regions  of  the  brain.  Figure  3  shows  pressure  at  air  cranium 
(skull)  interface.  The  blast  pressure  increases  about  3.6  times 


Figure  3:  Blast  pressure  Figure  4:  Uderwash  effect 
at  air  skull  interface  of  the  helmet 


Figure  4  is  from  a  blast  simulation  of  a  head  with  helmet. 
The  1 .3  cm  gap  between  helmet  and  head  allows  the  blast  wave 
to  wash  under  the  helmet.  When  this  “underwash”  occurs, 
geometric  focusing  of  the  blast  wave  causes  the  pressures  under 
helmet  cranium  subspace  to  exceed  those  outside  the  helmet  . 
This  is  in  turn  causes  additional  loads  on  the  top  portion  of  the 
skull  where  the  “underwash”  effect  is  most  dominant. 
Underwash  effect  is  however  not  dominant  in  front  regions  of 
the  skull.  Figure  5  (a)  and  (b)  compares  pressure  in  top  and 
front  regions  of  helmet  cranium  subspace.  The  pressure  in 
helmet  cranium  subspace  is  transferred  to  skull  and  brain. 
Figure  6  (a)  and  6  (b)  shows  pressure  in  top  and  front  regions  of 
the  brain  respectively.  As  seen  from  these  figures  helmet 
reduces  the  pressure  in  frontal  region  of  the  brain  however  it 
increases  the  pressure  in  top  regions  of  brain  due  to 
"underwash"  effect  occurring  in  top  regions  of  helmet  cranium 
subspace.  Similar  trend  is  observed  for  shear  stress  (mises)  in 
the  brain  and  contact  pressures  on  the  outer  surface  of  the  skull. 


(a)  7b) 

Figure  5:  Pressure  history  of  marked  locations  of  helmet 
cranium  subspace 


(a)  (b) 

Figure  6:  Average  pressure  in  the  brain,  (a)  Top  Region 
(b)  Front  Region 


CONCLUSIONS 

We  have  conducted  blast  simulations  to  understand 
response  of  head  with  protection  of  helmet.  Head  was  subjected 
to  0.3  MPa  blast  overpressure.  Our  simulation  shows  peak 
overpressure  rise  of  about  333.33  %  at  fluid  solid  interface.  We 
have  simulated  1.3  cm  gap  scenario  between  head  and  helmet. 
Helmet  protects  frontal  region  of  brain  as  it  attenuates  shock 
intensity  reaching  the  frontal  side  of  head.  On  the  other  hand, 
helmet  has  adverse  effect  in  protecting  top  regions  of  the  head 
due  to  ’underwash’  effect  in  top  regions  of  helmet  cranium  sub 
space. 
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